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Abstract We extend the theory of noise-induced phase synchronization to the case 
of a neural master equation describing the stochastic dynamics of an ensemble of un- 
coupled neuronal population oscillators with intrinsic and extrinsic noise. The master 
equation formulation of stochastic neurodynamics represents the state of each popula- 
tion by the number of currently active neurons, and the state transitions are chosen so 
that deterministic Wilson-Cowan rate equations are recovered in the mean-field limit. 
We apply phase reduction and averaging methods to a corresponding Langevin ap- 
proximation of the master equation in order to determine how intrinsic noise disrupts 
synchronization of the population oscillators driven by a common extrinsic noise 
source. We illustrate our analysis by considering one of the simplest networks known 
to generate limit cycle oscillations at the population level, namely, a pair of mu- 
tually coupled excitatory (E) and inhibitory (/) subpopulations. We show how the 
combination of intrinsic independent noise and extrinsic common noise can lead to 
clustering of the population oscillators due to the multiplicative nature of both noise 
sources under the Langevin approximation. Finally, we show how a similar analy- 
sis can be carried out for another simple population model that exhibits limit cycle 
oscillations in the deterministic limit, namely, a recurrent excitatory network with 
synaptic depression; inclusion of synaptic depression into the neural master equation 
now generates a stochastic hybrid system. 
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1 Introduction 

Synchronous oscillations are prevalent in many areas of the brain including sensory 
cortices, thalamus and hippocampus [1]. Recordings of population activity based on 
the electroencephalogram (EEG) or the local field potential (LFP) often exhibit strong 
peaks in the power spectrum at certain characteristic frequencies. For example, in the 
visual system of mammals, cortical oscillations in the y frequency band (20-70 Hz) 
are generated with a spatially distributed phase that is modulated by the nature of 
a visual stimulus. Stimulus-induced phase synchronization of different populations 
of neurons has been proposed as a potential solution to the binding problem, that 
is, how various components of a visual image are combined into a single coherently 
perceived object [2, 3]. An alternative suggestion is that such oscillations provide a 
mechanism for attentionally gating the flow of neural information [4, 5]. Neuronal os- 
cillations may be generated by intrinsic properties of single cells or may arise through 
excitatory and inhibitory synaptic interactions within a local population of cells. Ir- 
respective of the identity of the basic oscillating unit, synchronization can occur via 
mutual interactions between the oscillators or via entrainment to a common periodic 
stimulus in the absence of coupling. 

From a dynamical systems perspective, self-sustained oscillations in biological, 
physical and chemical systems are often described in terms of limit cycle oscilla- 
tors where the timing along each limit cycle is specified in terms of a single phase 
variable. The phase-reduction method can then be used to analyze synchronization 
of an ensemble of oscillators by approximating the high-dimensional limit cycle dy- 
namics as a closed system of equations for the corresponding phase variables [6, 7]. 
Although the phase-reduction method has traditionally been applied to deterministic 
limit cycle oscillators, there is growing interest in extending the method to take into 
account the effects of noise, in particular, the phenomenon of noise induced phase 
synchronization [8-15]. This concerns the counterintuitive idea that an ensemble of 
independent oscillators can be synchronized by a randomly fluctuating input applied 
globally to all of the oscillators. Evidence for such an effect has been found in exper- 
imental studies of oscillations in the olfactory bulb [11]. It is also suggested by the 
related phenomenon of spike-time reliability, in which the reproducibility of a single 
neuron's output spike train across trials is greatly enhanced by a fluctuating input 
when compared to a constant input [16, 17]. 

In this paper we extend the theory of noise-induced phase synchronization to the 
case of a neural master equation describing the stochastic dynamics of an ensem- 
ble of uncoupled neuronal population oscillators with intrinsic and extrinsic noise. 
The master equation formulation of stochastic neurodynamics represents the state 
of each population by the number of currently active neurons, and the state transi- 
tions are chosen such that deterministic Wilson-Cowan rate equations [18, 19] are 
recovered in an appropriate mean-field limit (where statistical correlations can be ne- 
glected) [20-23]. We will consider the particular version of the neural master equa- 
tion introduced by Bressloff [23], in which the state transition rates scale with the 
size N of each population in such a way that the Wilson-Cowan equations are ob- 
tained in the thermodynamic limit N — > oo. Thus, for large but finite N, the network 
operates in a regime characterized by Gaussian-like fluctuations about attracting so- 
lutions (metastable states) of the mean-field equations (at least away from critical 
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points), combined with rare transitions between different metastable states [24]. (In 
contrast, the master equation of Buice et at assumes that the network operates in a 
Poisson-like regime at the population level [21, 22].) The Gaussian-like statistics can 
be captured by a corresponding neural Langevin equation that is obtained by carry- 
ing out a Kramers-Moyal expansion of the master equation [25]. One motivation for 
the neural master equation is that it represents an intrinsic noise source at the net- 
work level arising from finite size effects. That is, a number of studies of fully or 
sparsely connected integrate-and-fire networks have shown that under certain condi- 
tions, even though individual neurons exhibit Poisson-like statistics, the neurons fire 
asynchronously so that the total population activity evolves according to a mean-field 
rate equation [26-30]. However, formally speaking, the asynchronous state only ex- 
ists in the thermodynamic limit N oo, so that fluctuations about the asynchronous 
state arise for finite N [31-34]. (Finite-size effects in IF networks have also been 
studied using linear response theory [35].) 

The structure of the paper is as follows. First, we introduce the basic master equa- 
tion formulation of neuronal population dynamics. We reduce the master equation to 
a corresponding neural Langevin equation and show that both intrinsic and extrinsic 
noise sources lead to multiplicative white noise terms in the Langevin equation. We 
then consider an ensemble of uncoupled neuronal populations each of which evolves 
according to a neural master equation. We assume that each population supports a 
stable limit cycle in the deterministic or mean-field limit. We apply stochastic phase 
reduction and averaging methods to the corresponding system of neural Langevin 
equations, following along similar lines to Nakao et al. [12], and use this to determine 
how independent intrinsic noise disrupts synchronization due to a common extrinsic 
noise source. (Previous studies have mostly been motivated by single neuronal oscil- 
lator models, in which both the independent and common noise sources are extrinsic 
to the oscillator. In contrast, we consider a stochastic population model in which the 
independent noise sources are due to finite size effects intrinsic to each oscillator.) 
We then apply our analysis to one of the simplest networks known to generate limit 
cycle oscillations at the population level, namely, a pair of mutually coupled exci- 
tatory (E) and inhibitory (/) subpopulations [36]. A number of modeling studies of 
stimulus-induced oscillations and synchrony in primary visual cortex have taken the 
basic oscillatory unit to be an E-I network operating in a limit cycle regime [37, 38]. 
The E-I network represents a cortical column, which can synchronize with other cor- 
tical columns either via long-range synaptic coupling or via a common external drive. 
In the case of an E-I network, we show how the combination of intrinsic independent 
noise and extrinsic common noise can lead to clustering of limit cycle oscillators due 
to the multiplicative nature of both noise sources under the Langevin approximation. 
(Clustering would not occur in the case of additive noise.) Finally, we show how a 
similar analysis can be carried out for another important neuronal population model 
that exhibits limit cycle oscillations in the deterministic limit, namely, an excitatory 
recurrent network with synaptic depression; such a network forms the basis of various 
studies of spontaneous synchronous oscillations in cortex [39—43]. We also highlight 
how the inclusion of synaptic depression into the master equation formulation leads 
to a novel example of a stochastic hybrid system [44]. 
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2 Neural Langevin equation 

Suppose that there exist M homogeneous neuronal subpopulations labeled i = 
1, ... ,M, each consisting of N neurons. 1 Assume that all neurons of a given sub- 
population are equivalent in the sense that the pairwise synaptic interaction between 
a neuron of subpopulation i and a neuron of subpopulation j only depends on i and 
j . Each neuron can be in either an active or quiescent state. Let Nf (t) denote the num- 
ber of active neurons at time t. The state or configuration of the full system (network 
of subpopulations) is then specified by the vector N(f) = (N\{t), Nz(t), . ■ ■ , A^m(0)> 
where each Ni(t) is treated as a discrete stochastic variable that evolves according 
to a one-step jump Markov process. Let P(n, t) — Prob[N(f) = n] denote the prob- 
ability that the full system has configuration n = (n\, n%, . . . , «m) at time t, t > 0, 
given some initial distribution P(n, 0). The probability distribution is taken to evolve 
according to a master equation of the form [20-23] 

dP(n, t) N N r 

*=lr=±l (1) 

-r t , r (n)PCn,0] 

with boundary condition P(n, t) = 0 if n, ■ — — 1 or n, ■ — N + 1 for some i. Here 
denotes the unit vector whose A:th component is equal to unity. The corresponding 
transition rates are chosen so that in the thermodynamic limit N — >• oo one recovers 
the deterministic Wilson-Cowan equations [18, 19] (see below): 

Tk-i (n) = a k n k , 

(2) 



T K+1 (n) = NF^Y, W kin,/N + I^j, 



where a k are rate constants, wu is the effective synaptic weight from the Zth to the 
kth population, and I\ are external inputs. The gain function F is taken to be the 
sigmoid function 

F(x)= — , (3) 

1 + e-y* 

with gain y and maximum firing rate Fq. (Any threshold can be absorbed into the 
external inputs /j.) Equation (1) preserves the normalization condition ^ Ml >o ■ ■■ 
Sn M >o P( n > 0 — 1 f°r all t > 0. The master equation given by equations (1) and (2) 
is a phenomenological representation of stochastic neurodynamics [20, 23]. It is mo- 
tivated by various studies of noisy spiking networks which show that under certain 



One could take the number of neurons in each sub-population to be different provided that they all scaled 
with N. For example, one could identify the system size parameter N with the mean number of synaptic 
connections into a neuron in a sparsely coupled network. 
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conditions, even though individual neurons exhibit Poisson-like statistics, the neu- 
rons fire asynchronously so that the population activity can be characterized by fluc- 
tuations around a mean rate evolving according to a deterministic mean-field equa- 
tion [26-29]. On the other hand, if population activity is itself Poisson-like, then it 
is more appropriate to consider an N -independent version of the master equation, in 
which NF — > F and yv/N w [21, 22]. The advantage of our choice of scaling 
from an analytical viewpoint is that one can treat N~ l as a small parameter and use 
perturbation methods such as the Kramers-Moyal expansion to derive a correspond- 
ing neural Langevin equation [45]. 

Multiplying both sides of the master equation (1) by n k followed by a summation 
over all configuration states leads to 



where the brackets (• • • > denote a time-dependent ensemble averaging over realiza- 
tion of the stochastic dynamics, that is, (/(n)> = £] n P(n, t) /(n) for any function 
of state f(n). We now impose the mean-field approximation (7fc r (n)) ^ 7]t jr ((n)), 
which is based on the assumption that statistical correlations can be neglected. Intro- 
ducing the mean activity variables x k — N~ l (n k ), we can write the resulting mean- 
field equations in the form 



Substituting for T kr using equation (2) yields the deterministic Wilson-Cowan equa- 
tions [19] 



Strictly speaking, the mean-field description is only valid in the thermodynamic limit 
N — >• oo, and provided that this limit is taken before the limit t — > oo [24]. In this 
paper we are interested in the effects of intrinsic noise fluctuations arising from the 
fact that each neural subpopulation is finite. 

Let us introduce the rescaled variables xi — n k /N and corresponding transition 
rates 




(4) 



r=±l 



-x k = N~ l [T k , + (Nx) - 7i_(iVi)] = H k (x). 



(5) 




(6) 




(7) 



Equation (1) can then be rewritten in the form 




(8) 



x P(x 



re k /N,t)-£2 k . r (x)P(x, t)]. 
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Treating x k as a continuous variable and Taylor expanding terms on the right-hand 
side to second order in N -1 leads to the Fokker-Planck equation 

<-| r>/ ^\ M „ 2 M <a2 

dP(X, t) 3 r T « X-v 3 r n 

^ = - E ^[ A *« p ( x < f) ] + yE ^2 [^(x)P(x, o] (9) 



with e — N~ 1/2 , 



and 



A*00 = Q/t,i(*)-n*-i(x) (10) 



i?t(i) = n tl i(x) + Q/ k -i(x). (ii) 

The solution to the Fokker-Planck equation (9) determines the probability den- 
sity function for a corresponding stochastic process X(r) with X(f) = . . . , 
•Xm(0)> which evolves according to a neural Langevin equation of the form 

dX k = A k (X) dt + 6y/B k (X) dW k (t). (12) 

Here W k (?) denotes an independent Wiener process such that 

(dW k (t)) = 0, (dW k (t)dWi(t)) = 8 kJ dt. (13) 

Equation (12) is the neural analog of the well known chemical Langevin equation [46, 
47]. (A rigorous analysis of the convergence of solutions of a chemical master equa- 
tion to solutions of the corresponding Langevin equation in the mean-field limit has 
been carried out by Kurtz [48].) It is important to note that the Langevin equation (12) 
takes the form of an Ito rather than Stratonovich stochastic differential equation 
(SDE). This distinction will be important in our subsequent analysis. 

The above neural Langevin equation approximates the effects of intrinsic noise 
fluctuations when the number N of neurons in each sub-population is large but finite. 
It is also possible to extend the neural Langevin equation to incorporate the effects of 
a common extrinsic noise source. In particular, suppose that the external drive I k to 
the A:th subpopulation can be decomposed into a deterministic part and a stochastic 
part according to 

2er Y k 

h=h k +—^m, d4) 



where h k is a constant input and £ (f ) is a white noise term, which is assumed to be 
common to all the neural subpopulations; the level of extrinsic noise is given by the 
dimensionless quantity a and ^Zj.=l Xk = !• Substituting for I k in equation (7) and 
assuming that a is sufficiently small, we can Taylor expand i to first order in a to 
give 

Q*,l(x) « f(J2 w ki*i + h k ^j + 2 ^f'{^ w km + h k ^(t). (15) 
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Carrying out a corresponding expansion of the drift function A k (x) then leads to the 
extended neural Langevin equation 

dX k = A k {X)dt + eb k (X)dW k (t)+aa k {X)dW{t), (16) 

where 

a k (x) - 2xk F' I ^2 Wkixi + hk I / y/~Fo, 

V I J (17) 

b k (x) = jB k (x), 

and dWit) — l-(t)dt is an additional independent Wiener process that is common 
to all subpopulations. We now have a combination of intrinsic noise terms that are 
treated in the sense of Ito, and an extrinsic noise term that is treated in the sense of 
Stratonovich. The latter is based on the physical assumption that external sources of 
noise have finite correlation times, so that we are considering the external noise to be 
the zero correlation time limit of a colored noise process. 



3 Stochastic synchronization of an ensemble of population oscillators 

In the deterministic limit (e, a — > 0) the neural Langevin equation (16) reduces to the 
mean-field equation (6). Suppose that the latter supports a stable limit cycle solution 
of the form x = x*(f) with x*(f + nT) = x*(t) for all integers t, where T is the period 
of the oscillations. The Langevin equation (16) then describes a noise-driven popula- 
tion oscillator. Now consider an ensemble of Af identical population oscillators each 
of which consists of M interacting sub-populations evolving according to a Langevin 
equation of the form (16). We ignore any coupling between different population os- 
cillators, but assume that all oscillators are driven by a common source of extrinsic 
noise. Introducing the ensemble label fj,, fi = 1, . , . , Af, we thus have the system of 
Langevin equations 

dxP> = AktX^A dt + eb k (xM) dW^\t) 

(18) 

+ (ja k (X (tJ - ) )dW(t), k=l,...,M. 

We associate an independent set of Wiener processes , k = 1 , . . . , M, with each 
population oscillator (independent noise) but take the extrinsic noise to be given by a 
single Wiener process W(t ) (common noise): 

(dW^\t)dW l (v) (t)) = S kJ S >l , v dt, (19) 

(dWjf-\t)dW(t)) = 0, (20) 

(dW(t)dW(t)) = dt. (21) 

Langevin equations of the form (18) have been the starting point for a number of re- 
cent studies of noise-induced synchronization of uncoupled limit cycle oscillators [9, 
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1 1-15]. The one major difference from our own work is that these studies have mostly 
been motivated by single neuron oscillator models, in which both the independent 
and common noise sources are extrinsic to the oscillator. In contrast, we consider a 
stochastic population model in which the independent noise sources are due to finite 
size effects intrinsic to each oscillator. The reduction of the neural master equation (1) 
to a corresponding Langevin equation (16) then leads to multiplicative rather than ad- 
ditive noise terms; this is true for both intrinsic and extrinsic noise sources. We will 
show that this has non-trivial consequences for the noise-induced synchronization of 
an ensemble of population oscillators. In order to proceed, we carry out a stochastic 
phase reduction of the full Langevin equations (18), following the approach of Nakao 
et al. [12] and Ly and Ermentrout [15]. We will only sketch the analysis here, since 
further details can be found in these references. We do highlight one subtle differ- 
ence, however, associated with the fact that the intrinsic noise terms are Ito rather 
than Stratonovich. 



3.1 Stochastic phase reduction 



Introduce the phase variable 6 e (— n, n] such that the dynamics of an individual limit 
cycle oscillator (in the absence of noise) reduces to the simple phase equation 9 -co, 
where co — 2jt/T is the natural frequency of the oscillator and x(f) = x*(9(t)). The 
phase reduction method [6, 7] exploits the observation that the notion of phase can 
be extended into a neighborhood M. C R 2 of each deterministic limit cycle, that is, 
there exists an isochronal mapping ^ : M. — >• [—it, n) with 9 — ^(x). This allows us 
to define a stochastic phase variable according to ©^(f) = <f(X^(f)) e [— it, jt) 
with X^(f) evolving according to equation (18). Since the phase reduction method 
requires the application of standard rules of calculus, it is first necessary to convert 
the intrinsic noise term in equation (18) to a Stratonovich form [25, 49]: 



4(xW)4i,(xW) 



GO 



dX 



dt 



2 

:bk(X W ) dW k W (t) + aa k (X w ) dW(t). 



(22) 



The phase reduction method then leads to the following Stratonovich Langevin equa- 
tions for the the stochastic phase variables ®w, /j. — 1, . . . , M [9, 12, 14]: 



Al 



k=l 



d®M = w + Z A .(0 (/i) ) fc J t(0 (/i) )3^(© ( ' /) ) dt 



(23) 



+ €b k (@ w ) dW k Ul \t) + aa k (e (p - ) ) dW 



Here Z k (6) is the Aith component of the infinitesimal phase resetting curve (PRC) 
defined as 



Z k (9) 



dx k 



(24) 



x=x*(6>) 
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with J2b=l Z k (9)A k (x*(9)) = a. All the terms multiplying Z k {9) are evaluated on 



the limit cycle so that 



a k (9) = a k (x*(9)), h(9) = b k (x*(9)), 

(25) 



db k (x) 

db k (8) 



dx k 



It can be shown that the PRC is the unique 2n -periodic solution of the adjoint linear 
equation [7] 

dZ M 

= -!>/*(** (0)2/(0. (26) 
1=1 

where A\ k — dAi/dx k , which is supplemented by the normalization condition 
J2 k Z k (t) dxt /dt = a>. The PRC can thus be evaluated numerically by solving the 
adjoint equation backwards in time. (This exploits the fact that all non-zero Floquet 
exponents of solutions to the adjoint equation are positive.) It is convenient to rewrite 
equation (23) in the more compact form 



d®^ = 



M 



dt 



(27) 



+ e /3/c(@ (M) ) dW^\t) + aa(@ (fi) ) dW(t), 



k=\ 



where 



M 



fi k m = Z k (9)b k (9), 0(0) = J2Zk(0)b k (6) db k (8), 

k=l 

M 

a{e)=Y J Z k {9)a k (9). 



(28) 



k=\ 



In order to simplify the analysis of noise-induced synchronization, we now convert 
equation (27) from a Stratonovich to an Ito system of Langevin equations: 

d& W =A w (G)dt + d$ w (®,t), (29) 

where {£ W(G, f)} are correlated Wiener processes with 0 = (0 (1) , . . . , ©^^.That 
is, 

M 

d^ w (<B, t) = 6^ft(0 w ) dW k W (t) + aa(& (fl) ) dW(t), (30) 

k=l 
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with (dfW(G.f)) = 0 and (d$ ( ^(0, t)d$ (v) (0, t)) = C (fiV) (®)dt, where 
C^ v) (0) is the equal-time correlation matrix 



M 



c^ v \e) = a\(d^)a{e^) + € 2 Y,h{6 w )h{B {v) )^v. (31) 



The drift term A^'(6) is given by 



1 

A w (0) = co- -S2{6 W ) + -B\e ( i A) ) (32) 



with 

B(9 (fl} ) = C M (e) 



M 

l2 , IST^rr, f n (iL\\l2 



■ o\a{e^)X + e 2 YXh(e w )Y 



(33) 



It follows that the ensemble is described by a multivariate Fokker-Planck equation of 
the form 

dt 3<9^ L J 

/i=i 

(34) 

M 



- y — — \c^ v \9)P(0,t)i 



2 ^ 3(9^3(9 

/x,y=l 

Equation (34) was previously derived by Nakao et al. [12] (see also [15]). Here, 
however, there is an additional contribution to the drift term A^' arising from the fact 
that the independent noise terms appearing in the full system of Langevin equations 
(18) are Ito rather than Stratonovich, reflecting the fact that they arise from finite size 
effects. 

3.2 Steady-state distribution for a pair of oscillators 

Having obtained the FP equation (34), we can now carry out the averaging procedure 
of Nakao et al, [12]. The basic idea is to introduce the slow phase variables ijr — 
f (Ar) ) according to O^ 1 = cot + f 1 and set Q(f, t) = P({cot + O^}, t). 
For sufficiently small e and a , Q is a slowly varying function of time so that we can 
average the Fokker-Planck equation for Q over one cycle of length T — 2jt/co. The 
averaged FP equation for Q is thus [12] 

dt 2 ^ 3ip 

M <35) 
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where 



and 



with 



£2(0) dd, 



£2= — / 

2tt J 0 

1 f n 
2tt J-* 



(36) 



(37) 



g(f) 



zjT •>-* k=i 



(38) 



)p k (o' + f)de'. 



Following Nakao et al. [12] and Ly and Ermentrout [15], we can now investigate the 
role of finite size effects on the noise-induced synchronization of population oscilla- 
tors by focussing on the phase difference between two oscillators. Setting Af — 2 in 
equation (35) gives 



d_Q 
dt 



= — £2 



32 32 



+ 



3^(1) 3^(2)] 2 



1 



9 



+ 



9^(1) 
d 2 



[a 2 g(0) + e 2 h(0)] 



Q 



9iA (1) 3V f(2) 
Performing the change of variables 

f = (fW+f ( V)/2, (p = f m - f m 
and writing Q{f (l) , iA (2) , 0 = *(V f > O*(0. 0 we obtain the pair of PDEs 

a* 6 2 _3* i r . , -.a 2 * 



and 



ao _ 3 
~dT~H 



^[a 2 (g(O)-g( ( p)) + e 2 h(Q)]0. 



These have the steady-state solution 
1 



To 



2;r cr 2 (g(0)- g((P)) + e 2 h(0) 



(39) 



where To is a normalization constant. 
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A number of general results regarding finite size effects immediately follow from 
the form of the steady-state distribution <t>o((/>) for the phase difference <f> of two pop- 
ulation oscillators. First, in the absence of a common extrinsic noise source (cr = 0) 
and e > 0, <3>o(<^>) is a uniform distribution, which means that the oscillators are com- 
pletely desynchronized. On the other hand, in the thermodynamic limit W->oowe 
have e = N~ 1 / 2 — ► 0 so that the independent noise source vanishes. The distribu- 
tion Oo(</>) then diverges at 6 — 0 while keeping positive since it can be shown that 
g(0) > g(@) [12]. Hence, the phase difference between any pair of oscillators accu- 
mulates at zero, resulting in complete noise-induced synchronization. For finite N, 
intrinsic noise broadens the distribution of phase differences. Taylor expanding g(<p) 
to second order in 0 shows that, in a neighbourhood of the maximum at (p — 0, we 
can approximate 4>o(</>) by the Cauchy distribution 

r' 

d>o(</>)~ 



(t> 2 <r 2 \g"(0)\/2 + h(0)/N 

for an appropriate normalization Y' Q . Thus the degree of broadening depends on the 
ratio 

^ _ MO) 

No 2 \g"{0)\ 

The second general result is that the functions a(6) and fik(&) that determine g(<f>) 
and h((p) according to equations (38) are nontrivial products of the phase resetting 
curves Zk(9) and terms atiO), bu(9) that depend on the transition rates of the original 
master equation, see equations (17), (25) and (28). This reflects the fact that both 
intrinsic and extrinsic noise sources in the full neural Langevin equation (18) are 
multiplicative rather than additive. As previously highlighted by Nakao et al. [12] 
for a Fitzhugh-Nagumo model of a single neuron oscillator, multiplicative noise can 
lead to additional peaks in the function g(4>), which can induce clustering behavior 
within an ensemble of noise-driven oscillators. In order to determine whether or not 
a similar phenomenon occurs in neural population models, it is necessary to consider 
specific examples. We will consider two canonical models of population oscillators, 
one based on interacting sub-populations of excitatory and inhibitory neurons and the 
other based on an excitatory network with synaptic depression. 

4 Excitatory-inhibitory (E-I) network 

4.1 Deterministic network 

Consider a network model consisting of an excitatory subpopulation interacting with 
an inhibitory subpopulation as shown in Figure 1(a). The associated mean-field equa- 
tion (6) for this so-called E-I network reduces to the pair of equations (dropping the 
bar on Jcj and setting M — 2) 



—— - -x E + F(w E exe + weixj + h E ), 
at 

dx[ 

—— - -xi + F(w IE x E + w u xi + hi), 
dt 



(40) 
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Fig. 1 Deterministic E-I network, (a) Schematic diagram of network architecture, (b) Phase diagram of 
two-population Wilson-Cowan model (40) for fixed set of weights wgg = 11.5, wjg = — wgj = 10, 
wji = —2. Also Fq — y = 1. The dots correspond to Takens-Bogdanov bifurcation points. 



where oiej — a — \ for simplicity. (We interpret a as a membrane time con- 
stant and take a~ l — 10 msec in physical units.) Also note that wee, w ie > 0 and 
wei, wii < 0. The bifurcation structure of the Wilson-Cowan model given by equa- 
tions (40) has been analyzed in detail elsewhere [36]. An equilibrium (x E ,Xj) is 
obtained as a solution of the pair of equations 

x| = F(w E ex e + weix* + hs), 

(41) 

x\ = F(wiex* e + wax*} + hi). 

Suppose that the gain function F is the simple sigmoid F(u) = (1 + e - ") -1 , that 
is, Fq—1 and y — 1 in equation (3). Using the fact that the sigmoid function then 
satisfies F' = F(l — F), the Jacobian obtained by linearizing about the fixed point 
takes the simple form 

j_ /-l + W E EX* E {\ -X* E ) WE1X E (\-X* E ) \ 

An equilibrium will be stable provided that the eigenvalues k± of J have negative 
real parts, where 

X±= i(TrJ±V[TrJ] 2 -4Detj). 

This leads to the stability conditions Tr J < 0 and Det J > 0. For a fixed weight ma- 
trix w, we can then construct bifurcation curves in the (x E , x*)-plane by imposing a 
constraint on the eigenvalues X± . For example, the constraint 

Tr J = -2 + w EE x* E {\ - x* E ) + wnx*(l - x*) = 0 

with Det J > 0 determines Hopf bifurcation curves where a pair of complex conjugate 
eigenvalues crosses the imaginary axis. Since the trace is a quadratic function of x%, 
x*, we obtain two Hopf branches. Similarly, the constraint Det J = 0 with Tr J < 
0 determines saddle-node or fold bifurcation curves where a single real eigenvalue 
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Fig. 2 Limit cycle in a deterministic E-I network with parameters wee = 11.5, wje = — WEI = 10, 
wji = —2, He = 0 and hj = —4. Also F(u) = 1/(1 +e - "). (a) Limit cycle in the (xg, xj) -plane, (b) Tra- 
jectories along the limit cycle for Xjj(r) (solid curve) and xj (f) (dashed curve). 

crosses zero. The saddle-node curves have to be determined numerically, since the 
determinant is a quartic function of x* E , x*. Finally, these bifurcation curves can be 
replotted in the (/z£, /z/)-plane by numerically solving the fix point equations (41) for 
fixed w. An example phase diagram is shown in Figure 1(b). 

We will assume that the deterministic E-I network operates in a parameter regime 
where the mean-field equations support a stable limit cycle. For concreteness, we 
take a point between the two Hopf curves in Figure 1 (b), namely, (h £,&/) = (0, —4) . 
A plot of the limit cycle is shown in Figure 2 and the components Zg, Z/ of the cor- 
responding phase resetting curve are shown in Figure 3. Note that both components 
are approximately sinusoidal so that the E-I network acts as a type II oscillator. 

4.2 Stochastic network and noise-induced synchronization 

Let us now consider an ensemble of uncoupled stochastic E-I networks evolving ac- 
cording to the system of Langevin equations (18) for M — 2 and k = E, I. (More 
precisely, each E-I network evolves according to a master equation of the form (1). 
However, we assume that N is sufficiently large so that the master equation can be 
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Fig. 3 Components Ze and Z; of phase resetting curve for an E-I network supporting limit cycle oscil- 
lations in the deterministic limit. Same network parameter values as Figure 2. 
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phase difference $ phase difference $ 

Fig. 4 Plot of periodic functions g and h for an E-I limit cycle oscillator with symmetric stochastic drive 
to excitatory and inhibitory populations (xe = XI = 1 /2). Same network parameters as Figure 2. 

approximated by the corresponding Langevin equation. This was also checked explic- 
itly in computer simulations.) Having numerically computed the phase resetting curve 
(Ze , Zj) and the solution on the limit cycle for the deterministic E-I network, we can 
then compute the functions g(<f>) and /;(</>) of the stationary phase distribution <t>o(0) 
according to equations (17), (25), (28) and (38). We plot these functions in Figure 4 
for the parameter values of the limit cycle shown in Figure 2, assuming symmetric 
common noise to excitatory and inhibitory populations. That is, xe — Xi — 1/2 in 
equation (17). It can be seen that the periodic function g is unimodal with g(0) > g(<f>) 
so that <£>o(<f>) is also unimodal with a peak at <f> — 0. 

The width and height of the peak depend directly on the relative strengths of the 
intrinsic noise e and extrinsic noise a . This is illustrated in Figure 5 where the ampli- 
tude a of the common signal is kept fixed but the system size TV is varied. Increasing 
N effectively increases the correlation of the inputs by reducing the uncorrected in- 
trinsic noise, which results in sharper peaks and stronger synchronization, see also 
Marella and Ermentrout [13]. We find that there is good agreement between our an- 
alytical calculations and numerical simulations of the phase-reduced Langevin equa- 
tions, as illustrated in Figure 6. We simulated the phase oscillators by using an Euler- 

0.26 „ , . . . . „ 




phase difference <)> 

Fig. 5 Probability distribution of the phase difference between a pair of E-I oscillators as a function of 
system size N for fixed extrinsic noise a = 0.08 with g, h given by Figure 4. Increasing N causes the 
curve to have a much sharper peak and much more synchronization. 
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Fig. 6 Probability distribution of the phase difference between a pair of E-I oscillators as a function 
of extrinsic noise strength a with g, h given by Figure 4. Solid curves are based on analytical calcula- 
tions, whereas black (gray) filled circles correspond to stochastic simulations of the phase-reduced (pla- 
nar) Langevin equations, (a) N = 10 5 , a = 0.01. The curve is very flat, showing little synchronization, 
(b) N = 10 5 , a = 0.08. Increasing a causes the curve to have a much sharper peak and much more syn- 
chronization. 



Mamyama scheme on the Ito Langevin equation (29). A large number A4 » O(l0 2 ) 
of oscillators were simulated up to a large time T (obtained by trial and error), by 
which time their pairwise phase differences had reached a steady state. As we were 
comparing pairwise phase differences each simulation gave us jM(Ai — 1) data 
points and we averaged over many simulations to obtain 10 6 data points for each dia- 
gram in Figure 6. These were then placed into 50 bins along [—it, it) and normalised. 
Also shown in Figure 6(b) are data points obtained from simulations of the full planar 
Langevin equations. Here computations were much slower so we only averaged over 
relatively few trials and thus the data is more noisy. Nevertheless a reasonable fit with 
the analytical distribution can still be seen. 

Nakao et al. have previously shown that in the case of Stuart-Landau or Fitzhugh- 
Nagumo limit cycle oscillators with both uncorrected and correlated extrinsic noise 
sources, parameter regimes can be found where the periodic function g has multi- 
ple peaks [12]. This can occur when higher harmonics of the phase resetting curve 
become dominant or when the common noise source is multiplicative. The presence 
of multiple peaks in g results in an ensemble of oscillators forming clustered states. 
Moreover, there are intermittent transitions between the clustered states induced by 
the uncorrected noise. In the case of stochastic E-I limit cycle oscillators, we were 
unable to find a parameter regime where g develops multiple peaks when the com- 
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Fig. 7 Periodic function g with multiple peaks when there is an asymmetry in the common extrin- 
sic noise source to the excitatory and inhibitory populations. Other network parameters are as in Fig- 
ure 2. (a) XE = 1/8. XI = 7/8 so that common stochastic drive is mainly to the inhibitory population, 
(b) XE = 7/8, // = 1/8 so that common stochastic drive is mainly to the excitatory population. 

mon extrinsic noise source is the same for both excitatory and inhibitory populations, 
that is, xe — XI — 1/2 in equations (14) and (17). However, multiple peaks can occur 
when there is an asymmetry between the excitatory and inhibitory stochastic drives, 
as illustrated in Figure 7. The corresponding stationary distribution $o(</>) for the 
phase differences <p also develops additional peaks, see Figure 8. When the common 
stochastic input is mainly presented to the inhibitory population, we find a peak at 
0 = 0 and smaller peaks at </> = ±2n /3. Consequently, the ensemble of oscillators 
tend to cluster in three regions around the limit cycle as shown in the inset of Fig- 
ure 8(a). On the other hand, when the stochastic drive is predominantly to the excita- 
tory population, we find a much sharper peak at 0 = 0 (compared to the symmetric 
case) and a small peak at <f> — it . However, the latter does not contribute significantly 
to the dynamics, so that the oscillators are strongly synchronized. 



5 Excitatory network with synaptic depression 

So far we have applied the stochastic phase reduction method to a two-population 
model consisting of mutually interacting excitatory and inhibitory populations. This 
E-I network is one of the simplest population models known to exhibit limit cycle 
oscillations in the deterministic limit, and forms the basic module in various studies 
of stimulus-induced oscillations and synchronization in visual cortex [37, 38]. An 
even simpler population model known to exhibit limit cycle oscillations is a recurrent 
excitatory network with synaptic depression. For example, Tabak et al. [39, 40] have 
analyzed Wilson-Cowan mean-field equations representing a recurrent excitatory net- 
work with both slow and fast forms of synaptic depression, and used this to model the 
dynamics of synchronized population bursts in developing chick spinal cord. These 
burst oscillations are more robust in the presence of an extrinsic noise source or some 
form of spatial heterogeneity within the network [50, 51]. An excitatory network with 
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Fig. 8 Probability distribution of the phase difference between a pair of E-I oscillators when there is 
an asymmetry in the common extrinsic noise source to the excitatory and inhibitory populations. Here 
N = 10 5 , a = 0.01 and other network parameters are as in Figure 2. Solid curves are based on analyt- 
ical calculations, whereas black filled circles correspond to stochastic simulations of the phase-reduced 
Langevin equations, (a) XE = l/8> XI = 7/8 so that common stochastic drive is mainly to the inhibitory 
population, (b) XE = 7/8, // = 1/8 so that common stochastic drive is mainly to the excitatory popula- 
tion. The insets show instantaneous distributions of the oscillators on the limit cycle. 



synaptic depression and extrinsic noise has also been used to model transitions be- 
tween cortical Up and Down states [41-43]. Here we will show how our analysis of 
noise-induced synchronization of population oscillators based on a Langevin approx- 
imation of a neural master equation can be extended to take into account the effects of 
synaptic depression. In addition to the relevance of synaptic depression in the gener- 
ation of neural oscillations, it is interesting from a mathematical perspective since the 
resulting master equation provides a novel example of a so-called stochastic hybrid 
system [44, 52]. 
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5.1 Deterministic network 



The mean-field equations for a homogeneous network with synaptic depression are 
taken to be of the form [39, 43, 53, 54] 

dx 

— = —x + F(qx + h), 

f < 42 > 
dq 

— = k + (\ — q) — k_xq, 
dt 

where we have set the membrane rate constant a = 1. The depression variable q(t) 
can be interpreted as a measure of available presynaptic resources at the population 
level, which are depleted at a rate k-x(t), which is proportional to the mean pop- 
ulation activity x(f), and are recovered at a rate k + . A fixed point (x*,q*) of the 
mean-field equation satisfies q* = k + /(k + + k-x*) with x* given by 

= Fl k + X * 



+ k-x* 

Suppose that the network operates in a parameter regime where there exists a unique 
fixed point. By linearizing about the fixed point and calculating the eigenvalues of 
the Jacobian, we can find regimes where the fixed point destabilizes via a Hopf bi- 
furcation leading to the formation of a limit cycle. An example bifurcation diagram 
is shown in Figure 9 with the depletion rate k- treated as the bifurcation parameter. 
Also shown is an example of a limit cycle for a particular value of k- . Given a limit 
cycle solution (x*(6), q*(0)) and associated isochronal function ^(x, q), we can nu- 
merically calculate the components of the corresponding phase resetting curve, which 




Fig. 9 Deterministic excitatory network with synaptic depression, (a) Bifurcation diagram showing solu- 
tions as a function of the depletion rate k— . Stable fixed point (thin solid curve) undergoes a Hop bifurca- 
tion at the points H, resulting in an unstable fixed point (dotted curve) and a stable limit cycle (thick solid 
curve), (b) Trajectories along the limit cycle for x (t) (solid curve) and q (t ) (grey curve). Parameter values 
are k+ = 0.02, y = 20, F 0 = 1 and h = -0.15 with k- = 0.1 in (b). 
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phase 9 

Fig. 10 Components Z x and Zq of phase resetting curve for an excitatory network with synaptic depres- 
sion. Same parameter values as Figure 9(b). 



are defined according to 

Z x (0) = 

Za(6) = 



dV(x,q) 



dx 
dV(x,q) 



x=x*(0),q=q*(6) 

(43) 



--x*(e),q=q*(8) 



dq 

together with the normalization condition 

co = Z x (G)A(x*(6),q*(6)) 

+ Z q (9)[k + (l-q*m) - k_x*(9)q*(0)]. 

The components of the phase resetting curve for the limit cycle shown in Figure 9(b) 
are plotted in Figure 10. As in the case of the E-I network, the resulting population 
oscillator is type II. However, the PRC is no longer approximately sinusoidal. 

5.2 Stochastic network and noise-induced synchronization 

We will now construct a stochastic version of the above population model by assum- 
ing that the population activity evolves according to a birth-death master equation 
with transition rates that depend on the depression variable q. (Both the determin- 
istic and stochastic models make a strong simplication by assuming that synaptic 
depression, which occurs at individual synapses, can be represented in terms of a sin- 
gle scalar variable q.) 2 Let N(t) denote the number of excitatory neurons active at 
time t, with P (n, t) — Prob[Af(f) = n] evolving according to the master equation (1) 
with M — 1 : 

dP(n,t) 

= T-{n+\)P(n+\,t) 

dt 



In order to relate the population depression variable q to what is happening at individual synapses, we 
label individual neurons within an excitatory network by the index a = 1 , . . . , N and assume that the 
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+ T+(n-l,t)P(n-l,t) (44) 
- [T_(n) + T + (n,t)]P(n,t) 

with P(— 1, t) = 0. The transition rates are taken to be of the Wilson-Cowan form 

7Li(n)=an, T+(n,t) = NF(q(t)n/N + 1), (45) 

where / is an external input, and g(f) satisfies 

dq , v N(r) 

^=fc f (l- ? (0)-*-X(09(0. *(0 = -j^. (46) 

The master equation (44) is non-autonomous due to the dependence of the birth rate 
T+ on q(t), with the latter itself coupled to the associated jump Markov process via 
the depletion rate k-X(t). Thus equation (46) is only defined between jumps, during 
which q evolves deterministically. 

The system defined by equations (44)-(46) is an example of a so-called stochas- 
tic hybrid model based on a piecewise deterministic process. This type of model 
has recently been applied to genetic networks [55] and to excitable neuronal mem- 
branes [44, 52, 56]. In the latter case, the hybrid model provides a mathematical 
formulation of excitable membranes that incorporates the exact Markovian dynamics 
of single stochastic ion channels. Moreover, the limit theorems of Kurtz [48] can be 



neurons are globally coupled. Suppose that the firing rate r a of the oth neuron evolves according to 



drg 
ill 



d 1ab 
dt 



b=\ 



= fc+(l - q ab ) - k-q ab r b . 



Summing the second equation with respect to b and dividing through by N leads to the following equation 
for i 



dq a 
dt 



= k + (\-q a )-k-q a x, 



provided that the following mean-field approximation holds: 



--xq a . 



If we then average the first equation with respect to a and again impose the mean field approximation, we 
see that 



dx 
dt 



+ N- l J2 F (laX + h). 
a=l 



Finally, noting that q a (t) — > q(t) for sufficiently large t (after transients have disappeared), we recover 
equations (42). In constructing a stochastic version of the network, we will assume that the above mean- 
field approximation still holds even though the activity variables are now random. See [70] for a recent 
discussion of the validity of mean-field approximations in a stochastic network model with synaptic de- 
pression. 
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adapted to prove convergence of solutions of the hybrid model to solutions of a cor- 
responding Langevin approximation in the limit N — >• oo and finite time, where N is 
the number of ion channels within the membrane [44, 52]. 

In the case of our stochastic hybrid model of an excitatory network with synaptic 
depression, we can heuristically derive a Langevin approximation by first carrying out 
a Kramers-Moyal expansion of the master equation (44). That is, setting x — n/N and 
treating x as a continuous variable, we Taylor expand the master equation to second 
order in \/N to obtain the Fokker-Planck equation 

dP(x,t) 9 r , e 2 9 2 r , 

J =- — [A(x,q)P(x,t)] + —— i [B(x,q)P(x,t)] (47) 
at ox 2 dx L 

with € = N~ x l 2 , 



and 



A(x, q) — F(qx + I) — ax 
B(x, q) — F(qx + /) + olx. 



(48) 
(49) 



The solution to the Fokker-Planck equation (47) determines the probability density 
function for a corresponding stochastic process X(t), which evolves according to the 
Langevin equation 

dX = A(X,Q)dt + eb(X,Q)dWi(t) (50) 

with b(x,q) 2 — B(x,q), Wi(t) an independent Wiener process and, from equa- 
tion (46), 

dQ = [k+(l-Q)-k-XQ]dt. (51) 

Following along similar lines to the E-I network, we can also include an extrinsic 
noise source by decomposing the drive to the excitatory population as 

I = h + -^=Ht), (52) 

where h is a constant input and f (f ) is a white noise term of strength a . Substituting 
for / in equation (50) and assuming that a is sufficiently small, we can Taylor expand 
to first order in a to give 

dX = A(X, Q)dt+eb(X, Q)dW\{t) + aa(X, Q)dW(t), (53) 

where W is a second independent Wiener process and 

a(x,q) = F'(qx+h)//Fo. (54) 

Suppose that we have an ensemble of excitatory networks with synaptic depression 
labeled fi — 1, . . . , J\f, each of which supports a stable limit cycle (x* (6) , q* (6)) in 
the deterministic limit. Carrying out a stochastic phase reduction along similar lines 
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to that of the E-I network, we obtain the following system of Stratonovich Langevin 
equations: 



for ij, = 1 , . . . , TV. Here 



»-yQ(e M ) 



J? + e^(0 ( ' i) )^w 1 (M) (O 



■CTa(0 (/i) )t/TV(O 



S2(0)=Z;t(0)i(0)aZ>(0), 
a(e) = Z x (9)a(6), 



P(9) = Z x (0)b(0), 



(55) 



(56) 



with 



a(6) =a(x*(0), q*(0)), b(0) = b(x*(0), q*(0)), 
db(x, q) 

x=x*(e),q=q*(9) 



db(&) 



fix 



(57) 



Thus, equations (55) and (56) correspond to equations (27) and (28) for M — 1 and 
Z\(6) = Z X (G). However, the functions Q(6), a (8), f}(6) implicitly depend on the 
dynamics of the depression variable as seen in equation (57). We can now write down 
the associated Fokker-Planck equation (34) and carry out the averaging procedure of 
Nakao etal. [12]. The final result of this analysis is the steady state distribution <t>o(</>) 
for the phase difference <p of any pair of oscillators given by equation (39) with 



km 



i r n 

- / a(e')a(6' + f)dd', 

TT J-„ 



(58) 
(59) 



and a, ft given by equation (56). An example plot of the periodic functions g(i/r), 
/z(i/f) for an excitatory network with synaptic depression is given in Figure 11. In 
Figure 12 we plot an example of the distribution <t>o illustrating how, as in the case 
of an E-I network, the synchronizing effects of a common extrinsic noise source are 
counteracted by the uncorrelated intrinsic noise arising from finite-size effects. 



6 Discussion 

In this paper we extended the theory of noise-induced synchronization to a stochastic 
Wilson-Cowan model of neural population dynamics formulated as a neural master 
equation. We considered two canonical network structures that are known to exhibit 
limit cycle oscillations in the deterministic limit; an E-I network of mutually inter- 
acting excitatory and inhibitory populations, and an excitatory network with synaptic 
depression. In both cases, we used phase reduction methods and averaging theory to 
explore the effects of intrinsic noise on the synchronization of uncoupled limit cycle 
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Fig. 11 Plot of periodic functions g and h for an excitatory network with synaptic depression. Same 
parameters as Figure 9(b). 



oscillators driven by a common extrinsic noise source. We achieved this by first ap- 
proximating the neural master equation by a corresponding neural Langevin equation. 
Such an approximation is reasonable for sufficiently large system size N, and pro- 
vided that there do not exist other stable attractors of the deterministic system [24]. 
One important consequence of intrinsic noise is that it broadens the distribution of 
phase differences. The degree of broadening depends on the term N~ l h(0), see equa- 
tion (39), where N is the system size and h(0) depends on the intrinsic dynamics of 
each uncoupled limit cycle oscillator. Another result our study is that the reduction of 
the master equation generates multiplicative rather than additive terms in the associ- 
ated Langevin equation for both intrinsic and extrinsic noise sources. Multiplicative 
noise can lead to clustering of limit cycle oscillators, as was demonstrated in the case 
of an ensemble of uncoupled E-I networks. 
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Fig. 12 Probability distribution of the phase difference between a pair of excitatory networks with synap- 
tic depression with g, h given by Figure 11. Solid curves are based on analytical calculations, whereas 
black filled circles correspond to stochastic simulations of the phase-reduced Langevin equation. Same 
parameters as Figure 9(b). 
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It is important to point out that the master equation formulation of stochastic neu- 
rodynamics developed here and elsewhere [21-24] is a phenomenological represen- 
tation of stochasticity at the population level. It is not derived from a detailed micro- 
scopic model of synaptically coupled spiking neurons, and it is not yet clear under 
what circumstances such a microscopic model would yield population activity consis- 
tent with the master equation approach. Nevertheless, if one views the Wilson-Cowan 
rate equations [18, 19] as an appropriate description of large-scale neural activity in 
the deterministic limit, it is reasonable to explore ways of adding noise to such equa- 
tions from a top-down perspective. One possibility is to consider a Langevin ver- 
sion of the Wilson-Cowan equations involving some form of extrinsic additive white 
noise [57, 58], whereas another is to view the Wilson-Cowan rate equations as the 
thermodynamic limit of an underlying master equation that describes the effects of 
intrinsic noise [20-23]. As we have highlighted in this paper, the latter leads to a 
multiplicative rather than additive form of noise. 

There are a number of possible extensions of this work. First, one could consider 
more complicated network architectures that generate limit cycle oscillations at the 
population level. One particularly interesting example is a competitive network con- 
sisting of two excitatory populations with synaptic depression (or some other form of 
slow adaptation) that mutually inhibit each other. Such a network has recently been 
used to model noise-induced switching during binocular rivalry [59-64]. Binocular 
rivalry concerns the phenomenon whereby perception switches back and forth be- 
tween different images presented to either eye [65, 66]. Experimentally, it has been 
found that the eye dominance time statistics may be fit to a gamma distribution, sug- 
gesting that binocular rivalry is driven by a stochastic process [67]. One possibility is 
that there is an extrinsic source of noise associated with the input stimuli. A number 
of recent models have examined dominance switching times due to additive noise in 
a competitive Wilson-Cowan model with additional slow adapting variables [61-63]. 
On the other hand, Laing and Chow [59] considered a deterministic spiking neuron 
model of binocular rivalry in which the statistics of the resulting dominance times 
appeared noisy due to the aperiodicity of the high-dimensional system's trajecto- 
ries. The latter is suggestive of an effective intrinsic noise source within a rate-based 
population model. A second extension of our work would be to introduce synaptic 
coupling between the limit cycle oscillators. For example, in the case of E-I networks 
such coupling could represent intracortical connections between columns in visual 
cortex [37, 38]. The effects of mutual coupling on noise-induced synchronization 
has been explored within the context of a pair of coupled conductance-based neu- 
rons [15]. Finally, the neural master equation has certain similarities to individual- 
based models in theoretical ecology, in particular, stochastic urn models of predator- 
prey systems [68, 69]. Given that predator-prey systems often exhibit limit cycle 
oscillations and receive extrinsic environmental signals, it would be interesting to 
extend our results on neuronal population oscillators to explore the effects of demo- 
graphic noise on the stimulus-induced synchronization of an ensemble of ecological 
communities. 
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